High resolution inverse synthetic aperture radar imaging of three-axis-stabilized space target by exploiting orbital and sparse priors
1. IntroductionWith the development of space technology, the number of space targets has increased dramatically, and the space environment is becoming more and more complex. In order to ensure the safety of space targets, inverse synthetic aperture radar (ISAR) imaging technology has been widely used for monitoring, tracking and identifying space targets in orbit.[1,2] The principle of two-dimensional (2D) ISAR imaging can be explained as the projection of the target onto the range–Doppler imaging plane. Usually, a high range resolution can be obtained by means of wideband transmitted signal, and a high cross-range resolution is acquired by using a large coherent integration rotation angle.[3,4] For high quality ISAR imaging, the translational motion must be compensated for accurately through range alignment[5,6] and phase adjustment.[7–10] However, we usually obtain the echoes of space targets in the case of low signal-to-noise ratio (SNR), and the range profiles are weak or completely submerged in noise. As a result, the performance of translational motion compensation with conventional methods degrades greatly. Fortunately, for a space target moving in orbit, its precise ephemeris can be used to compensate for the translational motion in practice. The precise ephemeris of a space target is calculated by using the telemetry data of the target, which is obtained by fusing the radars and optical sensors in a telemetry network. The residual error of precise ephemeris after subtracting an overall offset error could be able to reach a very small level and satisfy the requirement for range alignment. For phase correction, the required precision of ephemeris is higher than that of range alignment, which therefore must be considered in subsequent processing. Beyond that, it should be noted that in some imaging segments, the rotational motion of the target, which provides the cross-range resolution, is not a uniform rotation with a constant angular velocity, but the three-dimensional (3D) rotation with pitch, roll and yaw, which will cause blurring in the ISAR images. In this case, the nonlinear range migration and time-varying Doppler should be considered.
For a deliberately selected imaging segment which has a relatively stable imaging plane, the linear migration through resolution cells (MTRCs) caused by uniform planar rotational motion can be removed by the standard keystone transform algorithm (SKT),[11] and the sub echoes contain only a first-order phase item which can be processed with cross-FFT to obtain a well-focused ISAR image. Otherwise, for the imaging segments in which the spatial-variant imaging plane is not too severe, the space target can be treated as a general maneuvering target without any prior orbital information, and some maneuvering target imaging methods could be recommended to obtain an ISAR image. Most of the methods consider the sub echoes in a particular range cell as multicomponent linear frequency modulation (LFM) signals ignoring the slant-range MTRC caused by rotational motion. The time-frequency analysis[12–14] and multicomponent signal extraction techniques[15,16] have been introduced to achieve improved performances. For the purpose of reducing computational load, a modified keystone transform (MKT) is proposed for imaging the maneuvering targets,[17] which transforms the multicomponent LFM sub echoes into single-frequency signals, but the algorithm is usually efficient for the uniformly accelerated rotating target.
The recently developed sparsity-driven imaging is based on the fact that the target can be modeled to be sparse in the image domain. By using the sparse priors, the ISAR echoes of space target can be formulated as a sparse representation problem,[18–20] or an 
-norm optimization problem.[21–24] Furthermore, some novel sparsity-driven methods[25–28] have simultaneously taken into account the residual phase errors after translational motion compensation, but MTRC correction is not involved. In Ref. [29], a joint processing approach to range-invariant and range-variant phase errors is used for ISAR imaging with sparse aperture, but the approach may be applicable on the assumption that the rotational rate is uniform and the direction of effective rotation vector is constant, i.e., the pitch and roll motion of the target do not occur during the coherent processing interval (CPI).
If the imaging plane is severely spatial-variant, not only can the phase of sub-echoes not be approximated to a quadratic polynomial, but also the nonlinear MTRC both in slant- and cross-range dimensions must be solved. In this paper, we present a novel method for high resolution imaging of a three-axis-stabilized space target, which is more suitable for low SNR case by making full use of prior orbital information, and without making any assumptions for the equivalent rotational motion of space target during the CPI. The rest of this paper is organized as follows. In Section 2, the observation vector of radar to space target is obtained according to the orbital six elements. The change of scatterer coordinates is deduced through analyzing the three-dimensional rotation of the imaging plane, which is described by yaw angle, pitch angle and roll angle. In Section 3, the coarse translational motion compensation by using the precise ephemeris and the two-dimensional MTRC correction are proposed. In Section 4, the residual phase errors caused by the error of precise ephemeris are taken into consideration, and the sparsity-driven autofocus (SDA) method is adopted to adjust the residual phase errors and a high resolution ISAR image is achieved. Simulation results are shown in Section 5. Finally, some conclusions are drawn from the present study in Section 6.
2. The ISAR imaging geometry of three-axis-stabilized space targetThe position and attitude of an orbital target can be described precisely through two-line elements (TLE). In this section, we describe the scatterer positions at different observation times.
2.1. Representation of LOS in the ECEF coordinate systemThere are four coordinate systems commonly used for depicting the motion of space targets as shown in Fig. 1. Firstly we give the definitions of coordinate systems and the relationships between them.[30]
The coordinate system xyz is a geocentric observation coordinate system. Its origin is located in the center of the earth, x axis points to the space target from the geocentric point, the y axis is perpendicular to the x axis and points to the direction of target motion, obviously, the xy plane is the plane of the orbit, the z axis is normal to the xy plane.
The Earth-centered inertial (ECI) coordinate system XYZ is a non-rotating right-handed Cartesian coordinate system, the X axis points to the vernal equinox direction. The XY plane is the equatorial plane of earth, and the Z axis coincides with the earth axis of rotation and points northward.
The earth-centered earth-fixed (ECEF) coordinate system 
has a similar origin and direction of 
axis to the ECI system. The 
axis points to intersection of the prime meridian with the equatorial plane. The 
axis points to the intersection of the east longitude 90° with the equatorial plane. The coordinate system rotates together with the Earth.
The instantaneous imaging coordinate system 
is a right-handed coordinate system whose origin At is located in the mass center of the target, and subscript t is the observation time. The vt axis is the slant-range axis, pointing to the radar station C (fixed in ECEF), namely, the line-of-sight (LOS) direction. The ut axis is perpendicular to the vt axis in the imaging plane, and the wt-axis is orthogonal to the 
plane. The 
is an attitude-stabilized rotation unit vector of a three-axis-stabilized target. As the space target moves in orbit, the direction of 
changes over time.
The orbit of a space target can be determined by the six orbital elements which can be derived from a given TLE, and expressed as h(specific angular momentum), i (inclination), 
(right ascension (RA) of the ascending node), e (eccentricity), 
(argument of perigee) and θ (true anomaly). For a three-axis-stabilized target, its attitude remains unchanged with respect to the geocentric point, which means that the rotation angle caused by attitude stabilization of the target is equal to the true anomaly θ. According to the orbit equation,[31] at the observation time t, the coordinates of the mass center of target At in the geocentric observation coordinate system xyz can be determined by
where
h and
e are constants obtained from TLE, and
μ is the gravitational constant of the Earth, and its value is
.
According to the transformation relationship between the geocentric observation coordinate system and the ECEF system, we can obtain the coordinates of the mass center of target At in the ECEF system, which can be described by
where the rotation matrices
,
, and
express the rotation
ζ radians around the
x,
y, and
z axes respectively, and could be expressed as
Radar positions are usually given in the form of geographical longitude 
, latitude 
, and elevation 
, which can be transformed into the ECEF coordinate system by[32]
where
is the radius of the Earth, and its value is 6378.15 km generally.
Using Eqs. (3) and (5), we can obtain the radar’s unit observation vector relative to the target as
2.2. Scatterer position description in instantaneous imaging coordinate systemLet 
, 
, 
be the three axis unit vectors of coordinate system 
. As shown in Fig. 1, 
is the LOS vector, which is equal to 
. During imaging, the 
is time varying, and the velocity vector can be obtained by taking the slow time derivative of 
, which contains the tangential velocity vector and radial velocity vector. Thus the tangential velocity unit vector 
can be expressed as
The unit vector 
is normal to the 
plane, thus 
can be obtained by
As is well known, for a scatterer 
in Cartesian coordinate system uvw, the changes of its coordinates at time tk relative to initial position t0 can be described with a rotation transformation of coordinate system which contains three rotations around coordinate axes from initial time t0 to tk. Let ϕwk be the yaw angle rotating around the wk axis, which is normal to the imaging plane 
, ϕuk be the roll angle rotating around the uk axis, and ϕvk be the pitch angle rotating around the vk axis. We can obtain the rotation transformation matrix 
from Eq. (4) as follows:
 | |
For a three-axis-stabilized space target, the changes of axis vectors are caused not only by translational motion, but also by the attitude-stabilized rotation around vector 
, which is parallel to the normal vector of the orbital plane, going through the center of the target. The effects of attitude-stabilized rotation on the coordinate axis vectors can be eliminated with an anti-clockwise rotation transformation matrix 
, which has been discussed in Ref. [32]. Thus, the coordinate axis vectors at any observation time satisfy the following equation:
where
is the transformation matrix as follows:
where attitude-stabilized rotation vector
is known, and
are cosine and sine values of attitude-stabilized rotation angle
. Since the coordinate axis vectors
at any time can be determined by Eqs. (
6), (
7), and (
8), then we can calculate the rotation transformation:
The yaw angle, pitch angle and roll angle can be obtained as follows:
In the imaging coordinate system 
, the position change of scatterer p from time t0 to tk can be obtained by the following rotation transformation matrix:
where
and 
are the positions of scatterer
p at initial time
t0 and observation time
tk, respectively.
The vk axis points to the LOS direction, then the range projection value vpk of scatterer p onto the range axis can be obtained by Eqs. (9) and (14) as follows:
In Eq. (15), the product 
could be ignored, in which each factor is considered to be a small quantity. But the height item 
could cause the slant-range MTRC, according to the limitation 
. Assume that the radar transmits a signal with a bandwidth of 1 GHz, a scatterer with height 20 m in which there would not occur MTRC if the roll angle 
. During the CPI, the roll angle ξuk usually changes slightly with a level of 
, which will be verified in the following experiments with a real orbit. Thus the effect of height item could also be ignored, and equation (15) can be approximated as
3. Signal model of space target and MTRC correctionIt is assumed that the radar system transmits LFM signal. At time tk, for a scatterer 
in the target, the received signal after demodulation and range compression can be expressed as
where
, and
is the instantaneous range between the radar and the mass center of the target, which is the uniform part for each scatterer, and can be removed by translational motion compensation. However, in the case of low SNR situations, it is difficult to find some prominent points due to the effect of strong noise which would degrade the performance of translational motion compensation with conventional methods such as prominent points processing (PPP) and Doppler centroid tracking (DCT). Fortunately, for a space target in orbit, its precise ephemeris can be used to compensate for the echo phase partly, in other words, the precision of precise ephemeris could satisfy the requirement for envelope alignment, while its effect on echo phase cannot be ignored. Then, let
be the residual motion compensation error caused by the error in precise ephemeris, then the received signal after coarse translational motion compensation will be expressed as
Substituting Eq. (16) into Eq. (18), we can obtain the approximation:
where the
are re-represented as
for consistency. In Eq. (
19), the envelope term represents the position of the scatter, in which the coefficient of
vp is ignored, considering that the effect on
vp is less than half range cell. The phenomenon of the slant-range MTRC is mainly caused by the slant-range
up and the spatial-variant item
. The phase term provides the cross-range resolution, but also contains the phase errors, which is composed of slant-range varying phase errors, cross-range varying phase errors and spatial varying phase errors.
For an ideal imaging segment, where the pitch angle 
and roll angle 
are both very small, and the yaw angle 
is uniformly changed, the SKT can be used to correct the slant-range MTRC. In some segments which have a serious spatial-variant imaging plane, the pitch angle 
and roll angle 
cannot be ignored, and the yaw angle 
changes more complexly, the existing slant-range MTRC correction methods cannot be used directly. From Eq. (19), the slant-range MTRC can be described by using the three Euler angles which are obtained according to Eq. (13). Form this point of view, we can construct a relatively accurate scaling factor for the range-frequency to remove the slant-range MTRC. Firstly, equation (18) is transformed into the range-frequency domain with discrete form as
The scaling formula can be constructed according to Eq. (20), which is given by
where
PRT denotes the pulse repetition time (short as PRT),
is the average angular velocity during the CPI, and
denotes the nonuniform rotation factor. Substituting the scaling formula Eq. (
21) into Eq. (
20), we have
where
is the resampling slow-time. Here, the data
is transformed into
with a slow-time rescaled sampling, which can be realized by interpolation with Sinc or Knab kernels:
After applying inverse Fourier transform to Eq. (22) with variable f, the range profile without slant-range MTRC can be expressed as
where the second phase term on the right-hand side may cause the cross-range MTRC in ISAR image, and
represents the unknown spatial-variant phase error.
The ISAR image can be obtained by taking Fourier transform in cross-range direction after removing the cross-range MTRC from Eq. (24) by constructing the compensation item: 
, where 
, 
is the slant-range corresponding to operation unit, 
is the slant-range resolution cell, and 
is the unknown equivalent rotation center. Thus, the ISAR imaging can be rewritten in matrix form as follows:
where the operator
denotes Hadamard multiplication;
is a
matrix filled with the
K echoes, and
is
k-th echo data vector with a size of
;
is the
2D reflectivity matrix composed of
K1-D ISAR imaging
where the
znk is the reflectivity of the
-th cell,
and
respectively denote echo phase compensation matrix with a size of
and normalized Fourier transform matrix with a size of
, which are given by
and
where
.
The estimation of q can be obtained with a commonly accepted minimum entropy method, the objective function is expressed as
where
.
4. Sparsity-driven autofocus for high resolution ISAR imaging of a space targetAccording to the above-mentioned analysis, the residual phase errors can still cause defocus in cross-range direction even after previous processing. Let 
be the phase error of the kth echo, and 
be a 
diagonal matrix and represent the phase errors of K echoes, then the echoes with residual phase errors will be rewritten in the following form:
In order to improve the imaging resolution, the regularization-based imaging techniques can formulate the problem as an optimization problem which has been successfully applied to SAR and ISAR imaging. By exploiting the sparsity of 
, the sparsest solution of Eq. (27) can be converted into a sparsity-driven optimization as
where
is the
F-norm and denotes the fidelity term of the echoes,
is the
-norm and denotes the imposed sparsity,
λ is a regularization parameter and can be estimated using the proposed methods in Refs. [
26,
27], and [
29].
To avoid problems due to non-differentiability of the 
-norm at zero, a smooth approximation is used:[22]
where
ε is a nonnegative small constant, and the conjugate gradient algorithm (CGA) can be used to accelerate the convergence.
[23] Note that Eq. (
28) is a 2D optimization involving
F-norm, for finding the negative direction of the conjugate gradient vector of Eq. (
28), we need to convert it into a 1D problem that can be rewritten as
where
is the
vector obtained by stacking the columns of the echo matrix
S onto each other,
is the column-stacked version of the 2D reflectivity image
Z, and
Then, the conjugate gradient function of Eq. (30) is calculated through
Note that 
, then an iterative solver to Eq. (30) is presented through the following simplified formula:
where
and
g denotes the iteration number.
In fact, we do not need to construct the matrices 
and 
, equation (32) can be converted to 2D form as
where
is a
weight matrix formed by the diagonal elements of
.
For spatial-variant phase errors, the following cost function is minimized for each echo data vector 
:[25]
where
and
.
Minimizing the cost function (34) is equivalent to maximizing the item 
, in which the 
is a vector inner product, then the phase error estimation can be given in the closed forms of
and
In comparison with the phase error estimation employed in Ref. [27] and Ref. [29], the phase error estimation in Eq. (35) is concise, and the difference is the updating manner of 
. With the iteration format in Eqs. (33)–(35), we can obtain the high-quality inverse synthetic aperture radar image for estimating the phase errors. The iteration can stop whenever the relative change of the current estimation 
and previous estimation 
is smaller than a pre-determined threshold:
where
in this paper.
In order to illustrate the process clearly, the flow chart of the proposed imaging algorithm for a space target is given in the following Fig. 2.
5. Experimental resultsIn this section, we present experimental results to demonstrate the performance of the high resolution ISAR imaging method of a three-axis-stabilized space target. The radar works at a carrier frequency of 10.0 GHz and transmits LFM signal with a bandwidth of 400 MHz at a PRF of 50 Hz. Then the received signal is sampled at a rate of 500 MHz after mixing. The simulated target moves along a real orbit of Tiangong-1 satellite, and the echoes are generated according to the coordinate transformation model built in subSection 2.1. According to the previous analysis, the target attitude remains unchanged with respect to the geocentric point due to the attitude-stabilized motion of target, i.e., the range changes of scatterers, caused by its self-rotating motion are contained in phase history of echoes during the imaging time.
5.1. Spatial-variant property of imaging plane and MTRC correctionWe select the orbit of Tiangong-1 satellite as the simulation orbit, and its six orbital elements are provided by the Space Surveillance Network (SSN) of America in the form of a two-line element (TLE) as listed below.
1 37820U 11053A 16266.35688463 0.00025497 00000-024137-3 0 9991,
2 37820 042.7662 24.7762 0015742 351.0529 104.2087 15.66280400285808.
The epoch time of initial orbital elements is on 22 September 2016 at 08:33:54.8, the visible time period for radar placed in city Beijing is from 23:47:19 to 23:57:31; we chose three imaging segments which are shown in Fig. 3(a), where the red solid curve represents the orbit, and the blue lines denote the radar LOSs. The radar and space target are both in the Earth-fixed coordinate system. Three Euler angles of each segment are calculated according to the six orbital elements derived from the given TLE, which describe the spatial-variant property of the imaging plane in different imaging segments as shown in Fig. 3(b). With the same observation time (128 pulses), the segment 1 has a maximum yaw angle which is conducive to cross-range resolution and a minimum pitch angle which would not cause more serious MTRCs. By contrast, for constant observation time, the ISAR image reaches the best cross-range resolution in the imaging segment 1, the next is the resolution in segment 2, the worst is that in segment 3. Note that the roll angle (blue line in Fig. 3(b)) changes slightly in all segments, which implies that the approximation in Eq. (16) is reasonable.
The simulation model with 11 scatterers is shown in Fig. 4, and we generate 128 echoes for each imaging segment. The conventional range–Doppler (RD) imaging algorithm, which contains the range compression, the coarse phase correction, and the cross-range FFT, is adopted, and the imaging results without the MTRC correction of imaging segments 1–3, are shown in Figs. 4(a),4(d), and 4(g) respectively, from which we can see that the different degrees of the image distortion are caused by MTRCs: the more serious the spatial-variant property, the worse the imaging result is. Figures 4(b),4(e), and 4(h) show the imaging results of the three segments, obtained by using SKT to eliminate the MTRCs, and we cannot see the significant imaging quality improvement in these images. Figures 4(c),4(f), and 4(i) show the imaging results by using our proposed method in this paper to eliminate the MTRCs, and we can see that each image has a better focusing performance, which indicates that the MTRCs of different segments can all be eliminated effectively through our proposed method.
5.2. Elimination of residual spatial-variant phase errorsTo evaluate the performance of the phase error elimination, we further investigate the residual spatial-variant phase errors of an ideal scatterer in different imaging segments after MTRC correction. In order to obtain a constant cross-range resolution, two constant coherent integration rotation angles should be ensured in different imaging segments. In this simulation, two rotation angles of 1.685° and 2.815° are considered as shown in Table 1.
Table 1.
Table 1.
 Table 1.
Simulation orbital parameters of ISAR imaging.
.
| Segment |
1 |
2 |
3 |
| Start observation time |
23:50:15 |
23:53:15 |
23:56:15 |
| Pulse Number for integration angle of 1.685° |
188 |
460 |
1024 |
| Pulse Number for integration angle of 2.815° |
340 |
800 |
1860 |
| Table 1.
Simulation orbital parameters of ISAR imaging.
. |
The range compression and echo data interception are carried out first. The coarse translational motion compensation is then accomplished by utilizing the precise ephemeris derived from the instantaneous six orbital elements which can be provided by ground satellite station in reality. Finally, the sparsity-driven autofocus, detailed in Section 4, is used to eliminate the residual phase errors after slant- and cross-range MTRC corrections based on orbital priors. The blue lines in Figs. 5(a)–5(f) show the residual spatial-variant phase errors. Obviously they are those that are obtained by the longer radar observation, and the phase errors become more serious, which should not be ignored. After the phase errors are corrected by the proposed algorithm, the phase errors could be eliminated to a fairly low and satisfying level, no matter how big previous phase errors are, which are shown by the red lines in the following Fig. 5.
5.3. High resolution ISAR imaging of space target in different segmentsWe simulate a more complex space station with 307 scatterers shown in Fig. 6. The radar parameters and experiment orbital segments are the same as those listed in subSection 5.1. The observation time is chosen to ensure a constant cross-range resolution in each imaging segment, in this experiment the relative rotation angle is set to be 3.9°.
A complex Gaussian noise is added to the echo after range compression, and the SNR is set to be 5 dB. The conventional range alignment and coarse translational motion compensation are carried out first. Then, the SKT is utilized in slant-range MTRC correction, and the simulation results are shown in Figs. 7(a),7(d), and 7(g). Due to the spatial-variant property of imaging plane, the conventional method is not capable of eliminating the image distortion caused by MTRCs. The imaging results by using our proposed MTRC correction method are shown in Figs. 7(b),7(e), and 7(h). In contrast, we can see that the image has better focus performance. It is proved that the method of estimating the spatial-variant imaging plane of a three-axis-stabilized space target, using the target orbit information, is more effective in reducing the MTRCs, while the residual spatial-variant phase errors still blur the image to some extent. Finally, the sparsity-driven autofocus method is utilized to eliminate the residual phase errors, the corresponding imaging results are given in Figs. 7(c),7(f), and 7(i), with improved and satisfying quality.
5.4. Performance comparison among different autofocus imaging methods for space target under low SNRTo compare the sparsity-driven autofocus method analyzed in this paper with the existing autofocus techniques in low SNR situations, some simulation results are shown in this section. The imaging segment 3 shown in Fig. 3 is adopted, considering that the segment 3 has the most severely spatial-variant imaging plane. The echoes pre-processing method is the same as that described above. The SNR is set to be 0 dB, and a random phase error of uniform distribution in [
, 
] is added to sub echoes data. The imaging results are given in Figs. 8(a)–8(f) by using the dominant scatterer algorithm,[4] multiple prominent point processing (PPP),[4,8] Doppler centric tracking (DCT),[4,7] the weighted least-squares phase estimation (WLSPE),[3,10] phase gradient autofocus technique (PGA),[3,4,9] and the SDA, respectively. The comparison among the six images shows that the SDA method is capable of eliminating the spatial-variant error, and has a more satisfying imaging quality.
5.5. Measured data under low SNRA data set consisting of 400 pulses of the international space station (ISS) is recorded by an X-band ISAR experimental system, and the observation time is from 18:55:40 to 18:55:43 on 6 March 2010. After range compression on the raw echo data, the coarse translational motion compensation is applied to the data set, and the spatial-variant phase error is contained in each pulse. The aligned profiles are shown in Fig. 9(a). Due to the high radar cross section (RCS) of ISS, the signal-to-noise ratio of raw echoes reaches to 26 dB, which is high enough for proving available prominent scattering centers for conventional phase adjustment. To test the effectiveness of the proposed method under the low SNR, we add Gaussian complex noise into the echoes data to set the SNRs to be 15 dB and 5 dB, the polluted echoes are shown in Figs. 9(b) and 9(c). For comparison, we select the PGA to process the data set, which has the better performance in low signal-to-noise ratio case relatively, and generate the RD images shown in Figs. 10(a)–10(c). Then, we use the proposed method to obtain images as shown in the following Figs. 10(d)–10(f). It is clear that the proposed method can achieve a high quality ISAR image under the low SNR.
6. ConclusionsA practical high resolution ISAR imaging method for a three-axis-stabilized space target is proposed in this paper. Firstly, a set of spatial-variant angles of the imaging plane in instantaneous imaging coordinate system are calculated based on the orbital motion model of the three-axis-stabilized space target, with which the position histories of a scatterer in the target are accurately described. Then the slant- and cross-range MTRCs are corrected after the coarse compensation for translational motion. Finally, the sparsity-driven autofocus method is introduced to realize a fine compensation for residual phase errors and high-resolution ISAR image. The experimental results based on simulated and measured data illustrate that the proposed method is effective for the imaging segments in which the target moves with a three-dimensional rotation even under the condition of low SNR and the case of random phase errors.
